MODULE VEXCS_MOD
CONTAINS

SUBROUTINE VEXCS(KIDIA,KFDIA,KLON,KITT,KTILE,K_VMASS,LDINIT,&
 & PUMLEV,PVMLEV,PTMLEV,PQMLEV,PAPHMS,PGEOMLEV,PCPTGZLEV,&
 & PCPTS,PQSAM,PZ0MM,PZ0HM,PZ0QM,PZDL,PBUOM,PUCURR,PVCURR,PI10FGCV,&
 & YDCST,YDEXC,&
 & PCFM,PCFH,PCFQ,PKH,PKC)  

USE PARKIND1  , ONLY : JPIM, JPRB, JPRD
USE YOMHOOK   , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_THF   , ONLY : RVTMP2
USE YOS_EXCS  , ONLY : RCHBCD, DRITBL, RCHBBCD, RCHBB, RCHB23A, RITBL, &
 & RCHBA, RCHBHDL, RCDHALF, RCHETB, RCHBD, RCHETA, RCDHPI2, JPRITBL
USE YOS_CST   , ONLY : TCST
USE YOS_EXC   , ONLY : TEXC

! (C) Copyright 1990- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!     ------------------------------------------------------------------

!**   *VEXCS* - DETERMINES THE EXCHANGE COEFFICIENTS BETWEEN THE
!                 SURFACE AND THE LOWEST MODEL LEVEL WITH HELP OF
!                 STABILITY AS FUNCTION OF OBUKHOV-L.

!     Original A.C.M. BELJAARS       E.C.M.W.F.    26/03/90.
!     Modified A.C.M. BELJAARS   26/03/99 Tiling of the land surface.
!     Modified J. HAGUE          13/01/03 MASS Vector Functions       
!     Modified   P. Viterbo ECMWF 12/05/2005 Externalize SURF
!                                   (based on vdfexcs)
!     Modified A. Beljaars       30/10/2013 Change scaling of transfer coeff.
!     Modified A. Beljaars       02/02/2017 Introduction of tracer transfer coeff.
!     Modified M. Kelbling and S. Thober (UFZ) 11/6/2020 use of parameter values defined in namelist

!     PURPOSE
!     -------

!     DETERMINE EXCHANGE COEFFICIENTS BETWEEN THE SURFACE AND THE
!     LOWEST MODEL LEVEL

!     INTERFACE
!     ---------

!     *VEXCS* IS CALLED BY *SURFEXCDRIVER*

!     INPUT PARAMETERS (INTEGER):

!     *KIDIA*        START POINT
!     *KFDIA*        END POINT
!     *KLON*         NUMBER OF GRID POINTS PER PACKET
!     *KITT*         NUMBER OF ITERATIONS (BACK SUBST. IS ALWAYS DONE) [#]
!     *K_VMASS*      Controls the use of vector functions in the IBM scientific
!                     library. Set K_VMASS=0 to use standard functions

!     INPUT PARAMETER (LOGICAL):

!     *LDINIT*       IF .T. : ROUTINE GENERATES ITS OWN INITIAL GUESS  [#]
!                    IF .F. : ROUTINE STARTS WITH PZDL  AS INITIAL GUESS

!     INPUT PARAMETERS (REAL):

!     *PUMLEV*      X-VELOCITY COMPONENT AT T-1, lowest model level
!     *PVMLEV*      Y-VELOCITY COMPONENT AT T-1, lowest model level
!     *PTMLEV*      TEMPERATURE AT T-1, lowest model level
!     *PQMLEV*      SPECIFIC HUMUDITY AT T-1, lowest model level
!     *PAPHMS*      PRESSURE AT T-1, surface
!     *PGEOMLEV*    GEOPOTENTIAL AT T-1, lowest model level
!     *PCPTGZLEV*    DRY STATIC ENERGY, LOWEST MODEL LEVEL
!     *PCPTS*        DRY STATIC ENERGY AT THE SURFACE
!     *PQSAM*        SPECIFIC HUMIDITY AT THE SURFACE
!     *PZ0MM*        AERODYNAMIC ROUGHNESS LENGTH
!     *PZ0HM*        ROUGHNESS LENGTH FOR TEMPERATURE
!     *PZ0QM*        ROUGHNESS LENGTH FOR MOISTURE
!     *PZDL*         ZNLEV DEVIDED BY OBUKHOV LENGTH                 [#]
!     *PBUOM*        BUOYANCE FLUX AT THE SURFACE
!     *PUCURR*       OCEAN CURRENT X-COMPONENT
!     *PVCURR*       OCEAN CURRENT Y-COMPONENT
!     *PI10FGCV*     CONVECTIVE GUST VELOCITY

!     OUTPUT PARAMETERS (REAL):

!     *PCFM*         Scaled Transfer Coeff. momentum Rho*Cm*U      [kgm-2s-1]
!     *PCFH*         Scaled Transfer Coeff. heat     Rho*Ch*U      [kgm-2s-1] 
!     *PCFQ*         Scaled Transfer Coeff. moisture Rho*Cq*U      [kgm-2s-1] 
!     *PKH*          Scaled Transfer Coeff. heat         Ch*U      [m/s]
!     *PKC*          Scaled Transfer Coeff. tracer       Cc*U      [m/s]
!                    (for tracers, the coefficient describes the transfer 
!                     between the lowest model level and the roughness length
!                     for momentum, so it does not include the quasi-laminar 
!                     resistance)

!     REMARK: [#] UNUSED PARAMETERS IN TANGENT LINEAR AND ADJOINT VERSIONS
!     ------

!     METHOD
!     ------

!     THE ALGEBRAIC RELATION BETWEEN Z/L AND THE RICHARDSON NUMBER
!     IS SOLVED ITERATIVELY. THE STABILITY FUNCTIONS ARE THE SO-CALLED
!     PROFILE PSI-FUNCTIONS.
!     THE INITIAL GUESS (E.G. FROM PREVIOUS TIMESTEP) IS BACK-
!     SUBSTITUTED TO OBTAIN A SECOND APPROXIMATION. FURTHER ITERATION
!     IS DONE BY LINEAR INTER(EXTRA)POLATION (NEWTON'S  METHOD WITH
!     THE DERIVATIVE FROM SUCCESSIVE APPROXIMATIONS). IF NO INITIAL
!     GUESS IS PROVIDED, THE ROUTINE CAN PRODUCE ITS OWN. IN THE LATTER
!     CASE LDINIT=.T. HAS TO BE SPECIFIED AND IT IS RECOMMENDED TO
!     CHOOSE KITT=3. WITH INITIAL GUESSES FROM THE PREVIOUS TIME STEP
!     0 OR 1 ITERATION SHOULD BE SUFFICIENT.

!     ------------------------------------------------------------------

IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN)    :: KLON 
INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KFDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KITT 
INTEGER(KIND=JPIM),INTENT(IN)    :: KTILE 
INTEGER(KIND=JPIM),INTENT(IN)    :: K_VMASS
LOGICAL           ,INTENT(IN)    :: LDINIT 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PUMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PVMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PAPHMS(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PGEOMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCPTGZLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCPTS(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQSAM(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PZ0MM(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PZ0HM(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PZ0QM(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PZDL(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PBUOM(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PUCURR(:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PVCURR(:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PI10FGCV(:)
TYPE(TCST)        ,INTENT(IN)    :: YDCST
TYPE(TEXC)        ,INTENT(IN)    :: YDEXC
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCFM(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCFH(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCFQ(:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PKH(:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PKC(:) 
!*            LOCAL STORAGE
!             ----- -------
INTEGER(KIND=JPIM) :: JLEN

REAL(KIND=JPRB) ::  Z1DZ0M(KLON+MAX(K_VMASS,0))  ,Z1DZ0H(KLON+MAX(K_VMASS,0)),&
 & Z1DZ0Q(KLON+MAX(K_VMASS,0))  ,ZXLNM(KLON+MAX(K_VMASS,0)),&
 & Z1DZ0MD(KLON+MAX(K_VMASS,0)) ,Z1DZ0HD(KLON+MAX(K_VMASS,0)),Z1DZ0QD(KLON+MAX(K_VMASS,0)),&
 & ZXLNH(KLON+MAX(K_VMASS,0))   ,ZXLNQ(KLON+MAX(K_VMASS,0)),&
 & ZF2(KLON)             ,ZF3(KLON)             ,ZDU2(KLON+MAX(K_VMASS,0)),&
 & ZNLEV(KLON),&
 & Z1DZ1D(KLON+MAX(K_VMASS,0))  ,ZRICLS(KLON+MAX(K_VMASS,0)), ZETA2(KLON), ZETA3(KLON)
REAL(KIND=JPRD) :: ZDUMMY
REAL(KIND=JPRB) :: ZTMP1(KFDIA-KIDIA+1+MAX(K_VMASS,0))
REAL(KIND=JPRB) :: ZTMP2(KFDIA-KIDIA+1+MAX(K_VMASS,0))
REAL(KIND=JPRB) :: ZTMP3(KFDIA-KIDIA+1+MAX(K_VMASS,0))
REAL(KIND=JPRB) :: ZTMP4(KFDIA-KIDIA+1+MAX(K_VMASS,0))

INTEGER(KIND=JPIM) :: IRIB, JIT, JL

REAL(KIND=JPRB) :: ZA, ZAUX1, ZAUX2, ZB, ZCDNH,&
 & ZCDNM, ZCDNQ, ZCH, ZCM, ZCON1, ZCON2, ZCON3,&
 & ZCQ, ZDRORO, ZETA, ZETA1, ZHU, ZIPBL, &
 & ZL, ZPRH, ZPRH0, ZPRH1, ZPRM, ZPRM0, ZPRM1, &
 & ZPRQ, ZPRQ0, ZRIB1, ZRHO, ZUABS, ZWST2, &
 & ZX2, ZPRC0, ZPRC  

! for snow sublimation due to wind blowing
LOGICAL           :: LESNOWSUBL
REAL(KIND=JPRB)   :: ZSUBLSNW,ZUT,ZUV,ZASUB, ZBSUB
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!             INCLUDE STABILITY FUNCTIONS
!             ------- --------- ---------

#include "fcsvdfs.h"

!     ------------------------------------------------------------------

!*       1.   INITIALIZE CONSTANTS
!             ---------- ----------

IF (LHOOK) CALL DR_HOOK('VEXCS_MOD:VEXCS',0,ZHOOK_HANDLE)
ASSOCIATE(RD=>YDCST%RD, RETV=>YDCST%RETV, RG=>YDCST%RG, &
 & REPDU2=>YDEXC%REPDU2, RKAP=>YDEXC%RKAP, RPARZI=>YDEXC%RPARZI, &
 & RETACONV=>YDEXC%RETACONV )

ZCON1  =RVTMP2-RETV
ZCON2  =2.0_JPRB/3._JPRB
ZCON3  =RKAP**2

!     PBL HEIGHT FOR W* - EFFECT

ZIPBL=RPARZI

! Wind blowing snow sublimation is turned off by default
LESNOWSUBL=.FALSE. 
!     ------------------------------------------------------------------

!*       2.   COMPUTATION OF BASIC QUANTITIES
!             ----------- -- ----- ----------

IF(K_VMASS <= 0) THEN ! Not using Vector MASS

  DO JL=KIDIA,KFDIA

!             (W*)**2, WIND SHEAR, RICHARDSON NUMBER,

    IF (PBUOM(JL)  <=  0.0_JPRB) THEN
      ZWST2=0.0_JPRB
    ELSE
      ZWST2=(PBUOM(JL)*ZIPBL)**ZCON2
    ENDIF
    ZDU2(JL)=MAX(REPDU2,(PUMLEV(JL)-PUCURR(JL))**2&
     & +(PVMLEV(JL)-PVCURR(JL))**2+ZWST2) !+PI10FGCV(JL)**2)
    ZDRORO=2.0_JPRB*(PCPTGZLEV(JL)-PCPTS(JL))&
     & /(PCPTGZLEV(JL)+PCPTS(JL)-PGEOMLEV(JL))&
     & -ZCON1*(PQMLEV(JL)-MIN(PQSAM(JL),0.1_JPRB))
    ZRICLS(JL)=PGEOMLEV(JL)*ZDRORO/ZDU2(JL)

!             COMMON FACTORS IN NEUTRAL FORMULAE AND
!             DRAG COEFFICIENTS.

    ZNLEV (JL)=PGEOMLEV(JL)/RG+PZ0MM(JL)
    Z1DZ0M(JL)=ZNLEV(JL)/PZ0MM(JL)
    Z1DZ0H(JL)=ZNLEV(JL)/PZ0HM(JL)
    Z1DZ0Q(JL)=ZNLEV(JL)/PZ0QM(JL)
    Z1DZ1D(JL)=Z1DZ0M(JL)/(Z1DZ0M(JL)-1.0_JPRB)
    ZXLNM(JL) =LOG(Z1DZ0M(JL))
    ZXLNH(JL) =LOG(Z1DZ0H(JL))
    ZXLNQ(JL) =LOG(Z1DZ0Q(JL))
  ENDDO

ELSE ! Using Vector VMASS

  JLEN=KFDIA-KIDIA+K_VMASS-MOD(KFDIA-KIDIA,K_VMASS)
  DO JL=KIDIA,KFDIA
    ZTMP1(JL-KIDIA+1)=PCPTGZLEV(JL)+PCPTS(JL)-PGEOMLEV(JL)
    ZTMP3(JL-KIDIA+1)=ABS(PBUOM(JL)*ZIPBL)
  ENDDO 

  IF(KFDIA-KIDIA+1 /= JLEN) THEN
    ZTMP1(KFDIA-KIDIA+2:JLEN)=1.0_JPRB
    ZTMP3(KFDIA-KIDIA+2:JLEN)=1.0_JPRB
  ENDIF
  CALL VREC(ZTMP2,ZTMP1,JLEN)
  CALL VLOG(ZTMP4,ZTMP3,JLEN)
  DO JL=KIDIA,KFDIA
    ZTMP3(JL-KIDIA+1)=ZTMP4(JL-KIDIA+1)*ZCON2
  ENDDO
  CALL VEXP(ZTMP4,ZTMP3,JLEN)

  DO JL=KIDIA,KFDIA
!             (W*)**2, WIND SHEAR, RICHARDSON NUMBER,
    IF (PBUOM(JL)  <=  0.0_JPRB) THEN
      ZWST2=0.0_JPRB
    ELSE
      ZWST2=ZTMP4(JL-KIDIA+1)
    ENDIF
    ZDU2(JL)=MAX(REPDU2,(PUMLEV(JL)-PUCURR(JL))**2&
     & +(PVMLEV(JL)-PVCURR(JL))**2+ZWST2)
    ZDRORO=2.0_JPRB*(PCPTGZLEV(JL)-PCPTS(JL))*ZTMP2(JL-KIDIA+1)&
     &- ZCON1*(PQMLEV(JL)-MIN(PQSAM(JL),0.1_JPRB))
    ZTMP1(JL-KIDIA+1)=PGEOMLEV(JL)*ZDRORO

!             COMMON FACTORS IN NEUTRAL FORMULAE AND
!             DRAG COEFFICIENTS.

    ZNLEV (JL)=PGEOMLEV(JL)*(1.0_JPRB/RG)+PZ0MM(JL)
    Z1DZ0M(JL)=ZNLEV(JL)/PZ0MM(JL)
    Z1DZ0H(JL)=ZNLEV(JL)/PZ0HM(JL)
    Z1DZ0Q(JL)=ZNLEV(JL)/PZ0QM(JL)
    ZTMP2(JL-KIDIA+1)=Z1DZ0M(JL)-1.0_JPRB
  ENDDO

  IF(KFDIA-KIDIA+1 /= JLEN) THEN
    Z1DZ0M(KFDIA+1:JLEN)=1.0_JPRB
    Z1DZ0H(KFDIA+1:JLEN)=1.0_JPRB
    Z1DZ0Q(KFDIA+1:JLEN)=1.0_JPRB
    ZDU2  (KFDIA+1:JLEN)=1.0_JPRB
    ZDU2  (KFDIA+1:JLEN)=1.0_JPRB
    ZTMP2  (KFDIA-KIDIA+2:JLEN)=1.0_JPRB
  ENDIF
  CALL VDIV(ZRICLS(KIDIA),ZTMP1,ZDU2(KIDIA),JLEN)
  CALL VDIV(Z1DZ1D(KIDIA),Z1DZ0M(KIDIA),ZTMP2,JLEN)
  CALL VLOG(ZXLNM(KIDIA),Z1DZ0M(KIDIA),JLEN)  
  CALL VLOG(ZXLNH(KIDIA),Z1DZ0H(KIDIA),JLEN)
  CALL VLOG(ZXLNQ(KIDIA),Z1DZ0Q(KIDIA),JLEN)
ENDIF

!*       3. EXCHANGE COEFFIENTS DEPENDING UPON MONIN-OBUKHOV LENGTH
!           -------- ---------- --------- ---- ------------- ------

!           SOLVE Z/L ITERATIVELY
!           ----- --- -----------

IF(K_VMASS > 0) THEN
  CALL VREC(Z1DZ0MD(KIDIA),Z1DZ0M(KIDIA),JLEN)
  CALL VREC(Z1DZ0HD(KIDIA),Z1DZ0H(KIDIA),JLEN)
ELSE 
  DO JL=KIDIA,KFDIA
    Z1DZ0MD(JL)=1.0_JPRB/Z1DZ0M(JL)
    Z1DZ0HD(JL)=1.0_JPRB/Z1DZ0H(JL)
  ENDDO
ENDIF

!        3.1  FIRST AND SECOND GUESS TO INITIALIZE
!             ----- --- ------ ----- -- ----------

IF (LDINIT) THEN
  DO JL=KIDIA,KFDIA

!             FIRST GUESS: INTERPLOLATE ETA-TABLE FOR POS. RI
!             -----------

    IF (ZRICLS(JL)  >  0.0_JPRB) THEN
      ZRIB1=ZRICLS(JL)*ZXLNM(JL)/ZXLNH(JL)
      IRIB=INT(ZRIB1/DRITBL)+1
      IF (IRIB  >=  JPRITBL) THEN
        ZETA1=RITBL(JPRITBL)
      ELSE
        ZX2  = IRIB*DRITBL
        ZA   = (ZX2-ZRIB1)/DRITBL
        ZB   = 1.0_JPRB-ZA
        ZETA1= ZA*RITBL(IRIB) + ZB*RITBL(IRIB+1)
      ENDIF
    ELSE
      ZETA1=ZRICLS(JL)*ZXLNM(JL)/ZXLNH(JL)
    ENDIF

!             SECOND GUESS: WITH LOGARITHMIC TERMS
!             ------ -----

    ZETA2(JL)=ZETA1*ZXLNM(JL)*Z1DZ1D(JL)

  ENDDO
ELSE
  DO JL=KIDIA,KFDIA
    ZETA2(JL)=PZDL(JL)
  ENDDO
ENDIF

!        3.2  BACK SUBSTITUTION
!             ---- ------------

DO JL=KIDIA,KFDIA
  IF (ZETA2(JL)  >  0.0_JPRB) THEN
     ZDUMMY = ZETA2(JL)
     ZPRH1=PSIHS(ZDUMMY)
     ZDUMMY = ZETA2(JL)*Z1DZ0HD(JL)
     ZPRH0=PSIHS(ZDUMMY)
     ZDUMMY = ZETA2(JL)
     ZPRM1=PSIMS(ZDUMMY)
     ZDUMMY = ZETA2(JL)*Z1DZ0MD(JL)
     ZPRM0=PSIMS(ZDUMMY)
  ELSE
     ZDUMMY = ZETA2(JL)
     ZPRH1=PSIHU(ZDUMMY)
     ZDUMMY = ZETA2(JL)*Z1DZ0HD(JL)
     ZPRH0=PSIHU(ZDUMMY)
     ZDUMMY = ZETA2(JL)
     ZPRM1=PSIMU(ZDUMMY)
     ZDUMMY = ZETA2(JL)*Z1DZ0MD(JL)
     ZPRM0=PSIMU(ZDUMMY)
  ENDIF
  ZPRH=ZXLNH(JL)-ZPRH1+ZPRH0
  ZPRM=ZXLNM(JL)-ZPRM1+ZPRM0
  ZHU =ZPRM**2*Z1DZ1D(JL)/ZPRH
  ZETA3(JL) =ZRICLS(JL)*ZHU
  ZF2(JL)  =ZRICLS(JL)-ZETA2(JL)/ZHU
ENDDO

!        3.3  ITERATE KITT TIMES
!             ------- ---- -----

DO JIT=1, KITT
  DO JL=KIDIA,KFDIA
    IF (ZETA3(JL)  >  0.0_JPRB) THEN
       ZDUMMY = ZETA3(JL)
       ZPRH1=PSIHS(ZDUMMY)
       ZDUMMY = ZETA3(JL)*Z1DZ0HD(JL)
       ZPRH0=PSIHS(ZDUMMY)
       ZDUMMY = ZETA3(JL)
       ZPRM1=PSIMS(ZDUMMY)
       ZDUMMY = ZETA3(JL)*Z1DZ0MD(JL)
       ZPRM0=PSIMS(ZDUMMY)
    ELSE
       ZDUMMY = ZETA3(JL)
       ZPRH1=PSIHU(ZDUMMY)
       ZDUMMY = ZETA3(JL)*Z1DZ0HD(JL)
       ZPRH0=PSIHU(ZDUMMY)
       ZDUMMY = ZETA3(JL)
       ZPRM1=PSIMU(ZDUMMY)
       ZDUMMY = ZETA3(JL)*Z1DZ0MD(JL)
       ZPRM0=PSIMU(ZDUMMY)
    ENDIF
    ZPRH   =ZXLNH(JL)-ZPRH1+ZPRH0
    ZPRM   =ZXLNM(JL)-ZPRM1+ZPRM0
    ZF3(JL)=ZRICLS(JL)-ZETA3(JL)*ZPRH/(Z1DZ1D(JL)*ZPRM**2)
    IF (ABS(ZF3(JL)-ZF2(JL))  >  1.E-25_JPRB) THEN
      ZETA=ZETA2(JL)-ZF2(JL)*(ZETA3(JL)-ZETA2(JL))/(ZF3(JL)-ZF2(JL))
    ELSE
      ZETA=ZETA3(JL)
    ENDIF

!             SAFETY PROVISION FOR DIVERGING ITERATIONS

    IF (ZETA*ZETA3(JL)  <  0.0_JPRB) ZETA=ZETA3(JL)*RETACONV
    ZETA2(JL)=ZETA3(JL)
    ZF2(JL)  =ZF3(JL)
    ZETA3(JL)=ZETA
  ENDDO
ENDDO

!        3.4  COPY ZNLEV/L TO OUTPUT ARRAY
!             ---- ------- -- ------ -----

DO JL=KIDIA,KFDIA
  PZDL(JL)=ZETA3(JL)
ENDDO

!        4.   COMPUTE EXCHANGE COEFFICIENTS
!             ------- -------- ------------

DO JL=KIDIA,KFDIA

!             ASSUME MAXIMUM NOCTURNAL BOUNDARY LAYER HEIGHT
!             EQUAL TO RCHBHDL*L (5 TIMES OBUKHOV LENGTH)
!             AND UNIFORM PROFILES FROM H TO ZNLEV

  IF (ZETA3(JL)  >  RCHBHDL) THEN
    ZL=ZNLEV(JL)/ZETA3(JL)
    ZNLEV(JL)=RCHBHDL*ZL+MAX(PZ0MM(JL),PZ0QM(JL))
    ZETA3(JL)=ZNLEV(JL)/ZL
    Z1DZ0M(JL)=ZNLEV(JL)/PZ0MM(JL)
    Z1DZ0H(JL)=ZNLEV(JL)/PZ0HM(JL)
    Z1DZ0Q(JL)=ZNLEV(JL)/PZ0QM(JL)
    ZXLNM(JL) =LOG(Z1DZ0M(JL))
    ZXLNH(JL) =LOG(Z1DZ0H(JL))
    ZXLNQ(JL) =LOG(Z1DZ0Q(JL))
  ENDIF
ENDDO

IF(K_VMASS > 0) THEN
  CALL VREC(Z1DZ0MD(KIDIA),Z1DZ0M(KIDIA),JLEN)
  CALL VREC(Z1DZ0HD(KIDIA),Z1DZ0H(KIDIA),JLEN)
  CALL VREC(Z1DZ0QD(KIDIA),Z1DZ0Q(KIDIA),JLEN)
ELSE
  DO JL=KIDIA,KFDIA
    Z1DZ0MD(JL)=1.0_JPRB/Z1DZ0M(JL)
    Z1DZ0HD(JL)=1.0_JPRB/Z1DZ0H(JL)
    Z1DZ0QD(JL)=1.0_JPRB/Z1DZ0Q(JL)
  ENDDO
ENDIF

DO JL=KIDIA,KFDIA
  IF (ZETA3(JL)  >  0.0_JPRB) THEN
    ZDUMMY = ZETA3(JL)
    ZPRM1=PSIMS(ZDUMMY)
    ZDUMMY = ZETA3(JL)*Z1DZ0MD(JL)
    ZPRM0=PSIMS(ZDUMMY)
    ZDUMMY = ZETA3(JL)
    ZPRH1=PSIHS(ZDUMMY)
    ZDUMMY = ZETA3(JL)*Z1DZ0HD(JL)
    ZPRH0=PSIHS(ZDUMMY)
    ZDUMMY = ZETA3(JL)*Z1DZ0QD(JL)
    ZPRQ0=PSIHS(ZDUMMY)
    ZDUMMY = ZETA3(JL)*Z1DZ0MD(JL)
    ZPRC0=PSIHS(ZDUMMY)
  ELSE
    ZDUMMY = ZETA3(JL)
    ZPRM1=PSIMU(ZDUMMY)
    ZDUMMY = ZETA3(JL)*Z1DZ0MD(JL)
    ZPRM0=PSIMU(ZDUMMY)
    ZDUMMY = ZETA3(JL)
    ZPRH1=PSIHU(ZDUMMY)
    ZDUMMY = ZETA3(JL)*Z1DZ0HD(JL)
    ZPRH0=PSIHU(ZDUMMY)
    ZDUMMY = ZETA3(JL)*Z1DZ0QD(JL)
    ZPRQ0=PSIHU(ZDUMMY)
    ZDUMMY = ZETA3(JL)*Z1DZ0MD(JL)
    ZPRC0=PSIHU(ZDUMMY)
  ENDIF
  ZPRM   =ZXLNM(JL)-ZPRM1+ZPRM0
  ZPRH   =ZXLNH(JL)-ZPRH1+ZPRH0
  ZPRQ   =ZXLNQ(JL)-ZPRH1+ZPRQ0
  ZPRC   =ZXLNM(JL)-ZPRH1+ZPRC0
  ZRHO   =PAPHMS(JL)/(RD*PTMLEV(JL)*(1.0_JPRB+RETV*PQMLEV(JL)))  
  ZUABS  =SQRT(ZDU2(JL))
  ZAUX2  =ZUABS*ZRHO*ZCON3
  PCFM(JL)=ZAUX2/(ZPRM**2)
  PCFH(JL)=ZAUX2/(ZPRM*ZPRH)
  PCFQ(JL)=ZAUX2/(ZPRM*ZPRQ)
  PKH(JL) =PCFH(JL)/ZRHO
  PKC(JL) =ZUABS*ZCON3/(ZPRM*ZPRC)

! GA: snow sublimation due to wind blowing increases the transfer coefficient
!       for moisture in exposed snow and forest snow
        ! Also density always below 350
  IF (LESNOWSUBL) THEN
    IF (KTILE == 5 ) THEN
      ZUT=6.98_JPRB+0.0033_JPRB*(PTMLEV(JL)-245.88_JPRB)**2._JPRB
      ZUV=MIN(25._JPRB,SQRT(PUMLEV(JL)**2+PVMLEV(JL)**2)) ! Added a safety limiter to wind speed to avoid too large sublimation
      ZASUB=0.0018_JPRB
      ZBSUB=3.6_JPRB
      !*IF (ZUV>ZUT .AND. PRSN(JL) <= 350._JPRB)THEN
      IF (ZUV>ZUT)THEN
        ZSUBLSNW=ZASUB*(273.16/PTMLEV(JL))**4._JPRB*(ZUV/ZUT)**(ZBSUB-1_JPRB)
      ELSE
        ZSUBLSNW=0._JPRB
      ENDIF

      PCFQ(JL)=PCFQ(JL) + ZUABS*ZRHO*ZSUBLSNW ! rho*|U|*Cq + rho*|U|*Cq_snsub
    ENDIF
  ENDIF

ENDDO

END ASSOCIATE
IF (LHOOK) CALL DR_HOOK('VEXCS_MOD:VEXCS',1,ZHOOK_HANDLE)
END SUBROUTINE VEXCS

END MODULE VEXCS_MOD
